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This paper presents a unifying framework to uncertainty quantifica- 
tion for systems having polynomial response metrics that depend on both 
aleatory and epistemic uncertainties. The approach proposed, which is 
based on the Bernstein expansions of polynomials, enables bounding the 
range of moments and failure probabilities of response metrics as well as 
finding supersets of the extreme epistemic realizations where the limits of 
such ranges occur. These bounds and supersets, whose analytical struc- 
ture renders them free of approximation error, can be made arbitrarily 
tight with additional computational effort. Furthermore, this framework 
enables determining the importance of particular uncertain parameters ac- 
cording to the extent to which they affect the first two moments of response 
metrics and failure probabilities. This analysis enables determining the pa- 
rameters that should be considered uncertain as well as those that can be 
assumed to be constants without incurring significant error. The analytical 
nature of the approach eliminates the numerical error that characterizes 
the sampling-based techniques commonly used to propagate aleatory un- 
certainties as well as the possibility of under predicting the range of the 
statistic of interest that may result from searching for the best- and worst- 
case epistemic values via nonlinear optimization or sampling. 


I. Introduction 

Uncertainty Quantification (UQ) is the process of determining the effect of input uncer- 
tainties on response metrics of interest. Denote by p the input parameter vector whose value 
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is subject to uncertainty. These uncertainties may be classified as either aleatory , which are 
parameters subject to inherent and irreducible variability, or epistemic, which are reducible 
uncertainties resulting from a lack of knowledge. While the value assumed by aleatory vari- 
ables is expected to change for a given set of operating conditions (e.g., mass of fuel within 
a rocket), epistemic variables are unknown constants (e.g., mass of a structure). 

This paper studies the performance and reliability of a system whose response metrics are 
polynomial functions of the uncertain parameters. Regarding performance, the acceptability 
of the system depends upon the most likely outcome of a response metric. The performance 
analysis of a system consists of evaluating low order moments of the performance function 
for a given probabilistic model of the uncertainty. Regarding reliability, the acceptability of 
the system depends upon its ability to satisfy several design requirements simultaneously. 
These requirements, which are represented by a set of inequality constraints, depend on p. 
The system is deemed acceptable if all constraints are satisfied. The reliability analysis of a 
system consists of evaluating the probability of violating at least one of the requirements. 

A common approach to quantifying the effects of both aleatory and epistemic uncertain- 
ties is to perform nested sampling. This involves drawing samples of the epistemic variables 
in an outer loop and performing an UQ for the aleatory variables in an inner loop. In this 
fashion, ensembles of cumulative distribution functions, where each distribution represents 
the uncertainty generated by sampling over the aleatory variables, are generated. Nested iter- 
ation tends to be computationally expensive. Consequently, the nested sampling must often 
be under-resolved, particularly at the epistemic outer loop, resulting in an under-prediction 
of ranges of values of the response metric. Methods for propagating epistemic uncertainties, 
including interval-value probability, 2 second order probability, 3 imprecise probability theory 
and Dempster-Shafer theory of evidence, 10 suffer from such deficiencies. Methods that re- 
place the outer sampling loop with an optimization loop have been proposed in [1], Even 
though these methods are more efficient than sampling-based methods, their inability to 
guarantee convergence to the global optima may also result in under-estimates of the ranges 
of the statistics of interest. 

A key feature of this contribution is that the distinction between aleatory and epistemic 
variables implied by their definition is made consistently both qualitatively and quantita- 
tively. While aleatory uncertainties are manipulated according to long-standing concepts of 
probability theory, epistemic uncertainties are manipulated using properties of the Bernstein 
expansion of polynomials. In the proposed context each response metric becomes a random 
process. This family of random variables is parameterized by the value of the epistemic 
variables. In regard to the management of aleatory uncertainties, the framework proposed 
enables bounding tightly and rigorously the value of statistics supporting conventional UQ 
analyses (e.g. tight, formally verifiable ranges containing the failure probability range are 
calculated). In regard to the management of epistemic variables, the framework proposed en- 
ables bounding tightly and rigorously the epistemic realizations leading to the extreme values 
of such statistics (e.g., the realization of the epistemic variable leading to the largest failure 
probability). These bounds can be made arbitrarily small with additional computational 
effort. 

The strategies proposed yield exact results: the computed bound of a extreme statistic 
does not suffer from numerical or approximation error while the convergence to the criti- 
cal epistemic point realizing this extreme value is guaranteed. While standard probabilis- 
tic methods, such as polynomial chaos, Monte Carlo sampling, interval analysis, imprecise 
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probabilities, FORM, etc., cannot bound the approximation error present in their estimate 
of statistics, the identification of the extreme values attained by such statistics is out of their 
scope. As compared to methods based on interval analysis [5], whose results are also formally 
verifiable, the bounds proposed are better since they can always be made to converge to the 
set being bounded, e.g. bounds based on interval arithmetic suffer from irreducible conser- 
vatism when the requirement functions have repeated uncertain parameters. As compared 
to methods where the search for critical epistemic values is carried out using nonlinear pro- 
gramming or sampling, the proposed strategy eliminates the possibility of under-predicting 
the range of the statistic of interest, i.e., the search for the extrema of non-convex functions 
via nonlinear optimization may converge to non-global optima, thereby producing a UQ 
assessment that under-predicts the actual range. Furthermore, this framework enables the 
consideration of uncertainty models comprised of arbitrarily and possibly dependent aleatory 
variables, the substantial desensitization of the calculations from the assumed uncertainty 
model as well as the accommodation for changes in such a model with a small amount of 
computational effort. On the down side, as with all the UQ methods listed above, “the curse 
of dimensionality” restricts the applicability of the proposed strategies to systems with a 
moderate number of uncertainties. 

This paper is organized as follows. Basic concepts and notions are introduced in Sections 
II and III. Section IV presents strategies for UQ in the presence of aleatory uncertainties. 
This is followed by Section V where developments for UQ in the presence of epistemic and 
aleatory uncertainties are made. Low-dimensional, easily reproducible examples are used to 
benchmark the scope of framework. Finally, a few concluding remarks close the paper. 

II. Basic Concepts and Notions 

In a system that depends on the uncertain parameter p, response metrics are real- valued 
functions defined on a master domain PCM 3 . These metrics include performance functions 
and requirements functions. A performance function, such as cost, is a metric used to eval- 
uate a system’s performance. A reliability requirement is an admissible range of variation 
for a performance function, e.g., cost not exceeding the available budget. The functions 
prescribing reliability requirements will be called requirements functions. While a perfor- 
mance function will be denoted by g(p), where g : T> — >• M, the reliability requirements are 
prescribed as the vector inequality 3 " g(p) < 0, where g : T> — > ML 

The failure domain , denoted as T C M 3 , is comprised of the uncertain parameter realiza- 
tions that fail to satisfy at least one of the requirements. Specifically, the failure domain is 
given by 

S 

?= [J{P : 9i(p) > °} = {P ■ ™{P) > °}> (!) 

i= 1 

where w(p) = maxj<„{<7j(jp)} is the worst-case requirement function. The safe domain, given 
by S = U(J r ), where C(-) denotes the complement set operator given by C(Z) — V \ Z, 
consists of the parameter realizations satisfying all the design requirements. 

a Throughout this paper, it is assumed that vector inequalities hold component-wise, super-indices denote 
a particular vector or set, and sub-indices refer to vector components; e.g., p) is the zth component of the 
vector 
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Techniques for bounding T and S will be presented below. The resulting bounding sets 
are comprised of hyper- rectangles. The hyper-rectangle TZ G M s , whose “lower left” and 
“upper right” corners are at x and y, with x < y, is given by 

TZ(x, y) = {p : x < p < y} = (x u Vl ] x (* 2 , y 2 \ x • • • x (x s , y s ], (2) 

where the latter expression is the Cartesian product of intervals which exclude the left 
endpoint and include the right. Hyper-rectangles defined this way can be sub-divided into 
smaller hyper-rectangles without overlap, so that each point of the original hyper-rectangle 
falls into exactly one of the sub-dividing hyper-rectangles. Under these conditions, the 
larger hyper-rectangle is said to have been partitioned or subdivided into smaller hyper- 
rectangles. If p( TZ) = {7ii, ■ ■ ■ ,TZ t } is a pairwise disjoint collection of hyper-rectangles where 
li — Hi U . . . U lij, then p is a partition of TZ. A bisection-based subdivision of TZ C HU in 
the ith direction, i < s, is given by 

p(TZ(x, y)) = { TZ(x , x + w),K{y -w,y)}, (3) 

where w = [y 1 - aq, . . . , (y t - x t )/2, ...,y s -x s ]. 

The uncertainties considered in this paper are classified as either epistemic or aleatory. 
Denote by e G M e a subvector of p containing the epistemic variables, and by a G 1“ a 
subvector of p containing the aleatory variables, such that a + e = s. We will use non- 
probabilistic uncertainty models for e and probabilistic ones for a. Epistemic variables 
will be modeled by a bounded set. This set, called the support set, will be denoted as 
A e C M e . Each epistemic variable is modeled by providing an interval within which its value 
lies. Therefore, A e is the Cartesian product of these intervals. A probabilistic uncertainty 
model on the other hand, prescribes a measure of probability to each member of a set. 
Aleatory parameters will be modeled as random variables. This model is fully prescribed 
by the joint Probability Density Function (PDF) /(a) : A a C M a — y M, or equivalently, 
by the Cumulative Distribution Functions (CDF) F(a) : A a C M a — * [0,1]. Note that 
A a x A e C T> a xD f = T>. Not requiring the support set to fill up the master domain provides 
the freedom to change the uncertainty model without having to repeat calculations. 

The failure domain can also be described as 

r = U Me), (4) 

eSA e 

where 

F a (e) = {p= (a, e) : g(p) > 0}, (5) 

is a set value function of the epistemic uncertainty. 

The mean, variance and failure probability associated with the uncertainty model of the 
aleatory variables are 

E\g(p)\ = [ g(p) f (a)da, (6) 

J A a 

V[g{p)] = E[g\p )] - E[g{p)}\ (7) 

P[F] = [ f(a)da , (8) 

JTa 
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where E[-], V[-} and P[-] are the mean, variance and probability operators on the probability 
space A a . When some of the uncertain parameters are epistemic variables, these statistics 
become random processes parameterized by e. When all uncertain parameters are aleatory, 
these statistics are constants. The subscript “a” in T a and A a will be omitted hereafter 
since its meaning can be inferred from context. 

III. Bernstein Expansion 

The image of a hyper-rectangle when mapped by a multivariable polynomial is a bounded 
interval. By expanding that polynomial using a Bernstein basis over the rectangle, rigorous 
bounds to such an interval can be calculated by mere algebraic manipulations. Bernstein 
polynomials 6, 11 will be used for determining if a hyper-rectangle TZ is fully contained in 
the failure/safe domain or not. The outcome of the set containment test presented below 
depends exclusively on how much refinement of TZ the analyst is willing to perform. The 
refinement of TZ is determined by the number and size of subrectangles in a partition of TZ. 
Better refinements can always be used to render the set containment tests conclusive. The 
mathematical foundation of this approach is presented next. 

The Bernstein expansion first requires mapping 1Z into the unit hyper cube. Let KCD 
be an arbitrary hyper-rectangle in the master domain. Denote by u — U ( p ) an affine 
transformation that maps the hyper-rectangle TZ onto the unit hyper-cube U = TZ(0, 1). Let 
g be an arbitrary polynomial in M s . Then f{u) = giU^^u)) is a polynomial in U. Note 
that the extrema of g on TZ are identical to the extrema of / in U. 

For simplicity in the presentation we first consider a univariate polynomial. Since p, u 
and g are scalars here, we will represent them as p, u and g without the bold font. The 
transformation of g into the unit cube leads to 


f(u) = ^2 a i u \ 

i = 0 


whose Bernstein expansion is given by 


i = 0 

where 

B,» = ("W 

is the ith Bernstein polynomial of degree n (i.e., an element of the basis) and 


bi(TZ,g) = Y 


0 


=0 Vi/ 


" 3 % 


j = 0 


( 9 ) 


( 10 ) 


( 11 ) 


( 12 ) 


is the ith Bernstein coefficient. Some fundamental properties of this basis are ^ B™{u) = 1 
(normalization), 0 < Bf(u) < 1 for 0 < u < 1 (boundedness), and B^u) = (1 — u) > 

0 (symmetry). Some of the Bernstein coefficients assume the same value taken by the 
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polynomial at the vertices of a hyper-rectangle. This leads to the “free function evaluation 
property” , which is given by 

m = bo(K,g), (13) 

f(l) = b n (Jl,g). (14) 

The range enclosing property is described next. Suppose TZ is a hyper-rectangle and {bi(7Z, g ) : 
0 < i < n} are the Bernstein coefficients of g on TZ. The range enclosing property dictates 
that, for p e TZ, min o<i< n bi(7Z, g) < g(p ) < max 0 <i< n bi(TZ, g). Tighter bounds are obtained 
if TZ is subdivided. In particular, if p(JZ) = {TZ\, . . . , TZ t }, we have that, for all p e TZ, 

g(p, p) < g(p) < g(p, p% (15) 

where 

t 

g(p, p) = Y1 {Wh 9)} Kp ; Kj), ( 16 ) 

— L — ' 0 <i<n 

3 = 1 

t 

9(P,P) = ma x {bi{Kj,g)}I(p-,Tlj), (17) 

L ' 0<2<n 
3 = 1 ~ -■ 

where I(-;7Zj ) is the indicator function of TZp i.e., I(p;TZj ) = 1 if p G 7 Zj, I(p;7Zj) = 0, 
otherwise. Each of g and g is constant on each set 7 Zj. We call g and g the Bernstein lower 
and upper function bounds, respectively, of g. 

The multivariate polynomial case is considered next. Define the multi-index i to be a 
vector of non- negative integers of length s. The monomial ■ ■ ■ u\ s is abbreviated as 

u l . A s - variate polynomial can be represented as 

f(u ) = aiU *> ( 18 ) 

0<i<n 

where u G U. The Bernstein expansion of (18) is given by 

/(»)= E b i (K,g)B?(u), (19) 

0<i<n 

where 

BJ(u) = Bl'(u l )--Bl’(u,) (20) 

is the ith Bernstein polynomial of degree n and 

( p &) ^ 

bi(K,g) = J2 U.WB. V. (21) 

j<i \k= 1 \jj j 

for 0 < * < n, is the ith Bernstein coefficient. In this setting the free function evaluation 
property is 

/ ((ii/ni, . . . , i s /n s )) = 6 f ( 77, g), (22) 
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where i is an element of {0, rii} x • • • x {0, n p }. The range enclosing property is as follows: 
suppose 7 Z is a hyper-rectangle and {bi(JZ, g) : 0 < i < n} are the Bernstein coefficients 
of g on TZ. The range enclosing property dictates that, for p E TZ, min o<i<nbi(7Z, g) < 
g(p) < max 0 <i< n bi(JZ, g). As before, tighter bounds are obtained if Bernstein expansions 
over partitions of TZ are calculated. Specifically, 


for all p ETZ where 


g(p,p) < aip) < g(p,p), 


g(p, p) 

g{p , p) 


V min {bi{Kj,g)} I{p\Kj) r 

z — ■* 0 <t<n 
3 = 1 

t 

,' u h x i bi ( n j >g) } t (p ; n i)- ■ 

A ' 0 <t<n 
3 = 1 


(23) 

(24) 

(25) 


As before, g and g are constant on each sub-rectangle TZj. If p^ is the center of 7 Zj, 
these constant values are equal to g(p ^)) and ~g{p ( ^)- Since g is a continuous function on a 
connected set, we have 


Range (g) 


min g(p ) , max g (p) 
pgtz peiz 


min g(p ) , max g (p) 

P&7Z ~ p£TZ 


c 


(26) 


If the partition of 1Z is successively refined so that the maximum volume of the subsets 
approaches zero, the end points of the bounding interval converge to the global extrema of 

g- 

A strategy for bounding the values of p where the global extrema of a polynomial occur 
is presented next. Let p be a partition of 7 Z and let p * and p* denote the sets of global 
minima and global maxima, respectively, of g(p) over the hyper-rectangle. The lower and 
upper bounding functions in (24-25) can be used to calculate supersets of p * and p* . These 
supersets, denoted hereafter as V * and V*, are given by 

V * = U \ n 3 '■ g(P U) ) ^ min^(p)} D p *, (27) 

V * = U { n j : 9{P U) ) > g(p) j D p* . (28) 

3 ^ ' 


Therefore, the superset containing the point where g(p) attains its global minima, p *, is 
comprised of all rectangles where the minima of the upper bounding function g is larger 
than the lower bounding function g. Because the optimizations at the right hand side of the 
inequalities entail finding the extrema of piecewise constant functions over a finite partition 
of 7 Z, they can be determined exactly in a finite number of steps’ 3 . The supersets V * and 
V* approach p * and p* respectively as the partition of TZ becomes finer. The mathematical 

b Smaller supersets V * and V * result from replacing the arguments of the min and max operators in 
Equations (27) and (28) with a vector comprised of all free function evaluations, i.e., Equation (22) for all 
the elements of the partition. 
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foundation supporting the polynomial bounds guarantees that V * and V* will contain all 
global extrema regardless of their number and their location. Equations (26) yield lower 
and upper bounds of the extrema of g(p) while Equations (27) and (28) are supersets of the 
location of such extrema. 

The preceding analysis of a single polynomial function can be applied on a component by 
component basis to the vector g prescribing the failure domain. For a given partition of 7 Z, 
the preceding analysis is applied to each component of g to determine the bounding functions 
in (23). These functions can be used to calculate bounds of the worst-case requirement 
function in (1). The lower and upper bounding functions of w , denoted as w and w, are 
given by 


t 

™(p, p) = ^2 \ 

3 = 1 

( max g t (p b) ) 

\ l - 9 J 

| IfaKj), 

(29) 

t 

w(p, p) = 22 \ 

3 = 1 

fma xcnLp^) 

\ J 

'O'! 

.a, 

1 — i 

(30) 


where the indicator function / was dehned earlier. As before, w and w are piecewise constant 
on each member of the partition. The following theorems, which make use of these bounding 
functions, enable determining whether a set Q € is fully contained (or not) in the safe or 
failure domain. 

Theorem 1 (Set Containment in the Safe Domain). Let w(p ) be the worst-case require- 
ment function defined in (1) and pifLL) = {77i, . . . ,77. t } be a partition of the bounding set 'LL 
satisfying Q C 'LL. The set containment condition Q C S holds if 

max w(p, p) < 0. (31) 

p 

Furthermore, Q%S if there exists a k < v, a j < t, and a multi-index i e {0, rii} x • • • x 
{0, n p } such that 

bi{K j ,g k )> 0- (32) 

While Formula (31) results from using w(p) < w(p) for all p e 'LL in Equation (1), 
Formula (32) results from applying the free function evaluation property (22). 

Theorem 2 (Set Containment in the Failure Domain). Let w(p) be the worst-case require- 
ment function defined in (1) and pifLL) = {77i, . . . ,7 Z t } be a partition of the bounding set 'LL 
satisfying Q C 'LL. The set containment condition Q C T holds if 

min w(p, p) > 0. (33) 

p 

Furthermore, Q % LF if there exists a a k < v, a j < t, and a multi-index i e {0, rii} x • • • x 
{0,n s } such that 

hi (72-j i g if) < 0. (34) 

While Formula (33) results from using w(p) < w(p) for all p e 'LL in Equation (1), 
Formula (34) results from applying the free function evaluation property (22). 

The implicit formulation for calculating Bernstein coefficients proposed in [9] was adopted. 
This formulation is much more efficient than (12). Further efficiency can be realized by us- 
ing the subdividing logic in Equation (3) along with the algorithms of [6] that relate the 
Bernstein coefficients of a hyper-rectangle with those of its subsets. 
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IV. Uncertainty Quantification in the Presence of Aleatory 

Uncertainties 

In this section we develop strategies for bounding failure probabilities, means and vari- 
ances of response metrics that depend exclusively on aleatory uncertain parameters. 

A. Reliability Analysis 

By reliability analysis we refer to the bounding of the probability of failure. The key devel- 
opments in this section are the calculation of inner and outer bounding sets to the failure do- 
main and the calculation of their probabilities. These sets are comprised by hyper-rectangles 
whose membership into the bounding set is established using Theorems 1 and 2. Bernstein 
expansion-based techniques and interval analysis-based techniques 5 can be used to calcu- 
late these sets. The sets resulting from the latter technique however, are considerably more 
slack than those of the alternative technique for the same number of subsets. Furthermore, 
Bernstein bounds can always be made to converge to the function/set being bounded while 
bounds based on interval arithmetic may not, e.g. cases where the requirement functions 
have repeated parameters. 

This section presents an algorithm to generate and sequentially expand subsets of the 
failure and safe domains. These sets are unions of disjoint hyper-rectangles chosen from a 
partition p of T>. Let J 7 ™ 13 and 5 sub denote subsets (i.e., inner approximations) of the failure 
and safe domains formed from selected elements of p. Note that 0 C C T C C'(<S sub ) C 
T> and that the failure domain boundary, denoted as dF , lies in the region between the 
interiors of .F sub and d> sub . 

The sequences of inner bounding sets (<Sf b , d>f ub , . . .} and {.Ff ub , .Ff ub , . . .} are generated 
by the algorithm below. These sequences are made to converge to the domain being bounded. 
In particular, the algorithm iteratively generates indexed partitions Qi of T> and indexed 
sets (Sf ub , Jq sub and A* which are unions of hyper- rectangles from Qi, where <Sf ub is an inner 
approximation to the safe domain, J r f ub is an inner approximation to the failure domain, and 
A,; is a region comprised by the rectangles of Qi that are not in <Sf ub or Jq sub . Note that while 
Qi is a list of hyper-rectangles, S^ uh , J rt = uh and A* are sets comprised by the union of some 
of these rectangles. The algorithm proceeds by successively selecting each of the component 
hyper-rectangles TZ of A*. Then, Theorems 1 and 2 are used to determine if 7 Z is contained 
by the failure or safe domains. If the tests determine that TZ C S or TZ C J 7 , then TZ is 
removed from A* and added to <Sf ub or J 7 ™ 13 , respectively. If neither of these conditions is 
satisfied, 7 Z is subdivided, and the resulting sub-rectangles replace 7 Z in the partition. The 
algorithm terminates when the volume of A, : is sufficiently small. 

The algorithmic representation of this procedure is as follows. Let the inequality con- 
straint g(p) < 0 over p G V prescribe the system requirements. Set i — 1, Qi — {D}, 
Ai = T>, J 7 ™ 13 = 0, and iS™ b = 0. Pick a convergence criterion e > 0. 

1. Let L contain the elements of Q t comprising A*. 

2. Apply Theorems 1 and 2 to each hyper-rectangle in L without partitioning it, to 
determine which elements of L are contained in the safe domain and which ones are 
contained in the failure domain. Denote by U the list of elements contained by the 
safe domain, by V the list of elements contained by the failure domain, and by W the 
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list of elements that are not in U nor V. Furthermore, let U, V and W be the union 
of the elements in U, V and W respectively. 

3. Make Sf$ = Sf uh U U; = F| ub U V; and A i+1 = A ?: \ (U U V). 

4. If Vol(A,; + i) < e stop. Otherwise, make Q % +\ = (Qi \ W ) U p(W), increase i by one, 
and go to Step (1). 

As the number of iterations increases, <Sf ub and F) sub approach the safe and failure domain. 
Notice that the subdividing algorithm only partitions boxes whose containment in S or in T 
has not been established. Further notice that the closure of A, : not only covers the boundary 
of T but also approaches that boundary more and more closely as i increases. 

If a particular CDF defined over p e A C V is prescribed, bounds to the failure proba- 
bility, given by 

P[JT b ] < P[F\ < 1 - P[Si n % (35) 

can be readily calculated. These expressions are evaluated by adding up the probability 
of the rectangles comprising the bounding sets. Note that P[F| ub ] and P[5f ub ] are mono- 
tonically increasing functions of i; while P[A,] is a monotonically decreasing function of i. 
The choice of 77 made at step 1 of the algorithm implies that these bounds converge to 
Pf[P] as i increases. Therefore, all the conservatism in the bounds is reducible by additional 
computational effort. Further notice that once the bounding sets are available, failure proba- 
bility bounds corresponding to any distribution supported in a subset of V can be calculated 
efficiently. 

The values taken by w and w over the members of the partition p, which were required to 
compute the bounding sets iS sub and P sub , can be used to bound the CDF of the worst-case 
requirement function w. In particular, if p{V) = {Pi, . . . ,P*}, F w is bounded by 

Fw(x) = ^ F ^ x ) ^ Y. (36) 

jej( x) jej(x) 

where j(x) — {j : 1 < j < t,uJ(p ( - :j ' > ) < a:} and j(x) — {j : 1 < j < t,w(p^' ) ) < x} for all x 
in the range of iu. A few manipulations lead to 

1-P»(0)<P[P]<1-Fht(0). (37) 

This expression is equivalent to (35). The evaluation of (36) for each member of the partition 
sequence Qi resulting from the algorithm above yields a sequence of CDF bounds. These 
CDFs become tighter near w = 0 as i increases. 


Example 1: Consider the requirement functions 


9 1 

92 


2 4, 4 2 o 2 2 

P 1 P 2 + P 1 P 2 ~ 3piP 2 - PiP 2 + 


pIpI 

2 


- p{pl + 3£hJ>2 + 


Plgf 

10 


gi +pI 
200 
9 

~ 10 ’ 


7 

Too’ 


(38) 

(39) 


dehned over the master domain V G [—2, 2] x [—2, 2], The parameters p 1 and p 2 are aleatory 
uncertain parameters, i.e., a = {fq,^} an d e = 0, distributed according to the PDF 


fpiP2 (Pl>P 2 ) 


cos 2 {PiPQ 

8 + Si(8) ’ 


( 40 ) 
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where Si(-) is the sine integral. Note that the uncertainty models of p 1 and p 2 are strongly 
dependent. For instance, the conditional probability density function of p 1 given p 2 is 
uniform at p 2 = 0, and it is tri-modal at p 2 = —2. 

Figure 1 shows the set approximations J-f ub and S^ ub for a fixed value of i. Boxes 
comprising J r ? ub are colored in red, those comprising S^ ub are colored in green, and boxes 
comprising A* are colored in white. Note that A, : is a tight approximation of dJ-. Further 
notice that the density of boxes per unit of area increases with the closeness to dJ - . The 



Figure 1. 


F sub (red) and S sub (green) and failure domain boundary (thick line). 


corresponding failure probability bounds are 0.492 < P[J-} < 0.571. This probability cannot 
be calculated by standard sampling techniques. 

B. Performance Analysis 

Lower and upper bounds of moments of any performance function that depends on aleatory 
variables can be calculated from Equation (23). If the performance function g is a single 
polynomial in p and the aleatory parameters are independent random variables, moments 
of any order can be calculated exactly. This is not the case when either g is a continu- 
ous piecewise polynomial function, or when the aleatory variables are not independently 
distributed. 

In this section we develop techniques for bounding E[g(jp)\ and V[g(p)\ when g is piece- 
wise polynomial, and p is arbitrary distributed within a bounded set. The application of 
such techniques requires the calculation of the bounding functions in Equation (23). Even 
though the approach is applicable to any partition of the master domain, tighter bounds are 
obtained when the partition is finer in regions where either g — g or E[g — g\7Z] are large. 

In regard to the mean, the application of the expected value operator to Equation (23) 
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yields 


t t 

Y.sJP W )P[Kj\ < E[s] < J29(P U) )Pm- (41) 

3 = 1 J =1 

Bounds on the variance are considered next. Interval arithmetic will be used to simplify 
the notation. The starting point is the equation V[g\ = E[g 2 ] — E[g] 2 and the interval 
inclusion g G [g, g \ . We recall the formulae for interval subtraction and interval squaring: 

[a, b] — [c, d] = [a — d, b — c] and 


[ la 2 ,b 2 } 

M 2 = [& 2 ,a 2 ] 

[ [0, max {a 2 , b 2 }] 


if 0 < a, 
if b < 0, 

otherwise 


(42) 


In this context, g 2 G [g,g] 2 and E[g 2 ] G E [g,g] 2 = [/i, v}. The evaluation of this equation 


yields 

M = £ 2 (p U) ) 2 -P[R,] + J2 5(P U) ) 2 r[Kj], (43) 

jet/ jeV 

-= E 9(P U) ) 2 r[TCj] + E 9(P U) ) 2p [ K l]. ( 44 ) 

jeuuw jev ux 

where 

U = {1 < j < t : g(p^) > 0}, (45) 

V = {1 < j < t : g(p u) ) < 0}, (46) 

W = {1 < j < t : g(p u) )g(p u] ) < 0, \g(p {3) )\ < |^(p 0) )|)}, and (47) 

X = {1 < j <t :j £ (UUVUW)}. (48) 


The desired bounds are given by V[g] G E [[g,g] 2 ] — [E[g], [7y] ] 2 , whose evaluation yields 

max {0,o:} < V[g] < f3 , (49) 

where 

[ l-i ~ E[g} 2 if E[g\ > 0, 

ot — I g — E[g] 2 if E[g\ < 0, (50) 

[p — max {E{g] 2 , E\g] 2 } otherwise, 


v - E[g\ 2 if E\g\ > 0, 

P — | v — E[g] 2 if E[g\ < 0, (51) 

I v otherwise . 


Example 2: Given the same problem statement of Example 1, we want to bound E\w] 
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Figure 2. 


jrsub ( rec i) an d gsub (g reen ) a nd f a ilure domain boundary (thick line). 


and V[w \ , where w is the worst-case requirement function. Figures 2 and 3 show the set 
approximations J r f ub and Sf ub for a fixed value of i corresponding to two partition logics. 
The partition used to generate Figure 2 only subdivides boxes where g — g is in the top 50%. 
This partition leads to 3.21052 < E[w) < 4.71449 and 80.97456 < V[w\ < 145.11712. In this 
case, the density of boxes per unit of area increases with the separation between the upper 
and lower bounds of w. This separation is an indicator of the size of the range of w over a 
box. Consequently, the size of a box tends to be inversely proportional to the magnitude of 
the gradient of iv within the box. Note that the determination of whether a given box is in 
the failure or safe domain, therefore the coloring of the boxes in the figure, is inconsequential. 

The partition used to generate Figure 3 only subdivides boxes where E[g — g\7Z] is in the 
top 50%. This yields 3.40700 < E[w] < 4.50007 and 82.9050 < V[w] < 143.5641. This set of 
bounds are tighter than those from Figure 2 because they take into account the uncertainty 
model. In this case, the size of a box tends to be inversely proportional to the magnitude of 
the probability-weighted gradient of w within the box. 

V. Uncertainty Quantification in the Presence of Aleatory and 

Epistemic Uncertainties 

In this section we develop strategies for bounding failure probabilities, means and vari- 
ances of response metrics that depend on both aleatory and epistemic parameters. Recall 
that the propagation of probabilistic uncertainty from the aleatory variables to the response 
metrics yields a random process whose statistics are parametrized by the epistemic variables. 
Note that for a fixed value of the epistemic variable, the UQ methods of Section IV apply. 
While for the aleatory-only case the statistics take on a constant value, taking into account 
the epistemic intervals of uncertainty spreads each aleatory statistic over its own interval of 
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Figure 3. 


jTsah ( rec i) anc i gsub (g reen ) a nd f a ilure domain boundary (thick line). 


uncertainty. Each value in that interval is a possible realization of the statistic. The values 
of the epistemic realizations prescribing the interval’s limits will be referred to as best-case 
and worst-case epistemic values. These values will be denoted hereafter as e* and e*. In 
this section we use the Bernstein expansion approach to bound the range statistics of the 
response metric and calculate supersets of the corresponding e* and e*. The mathematical 
background for this is presented next. 

Denote by Proj a (A) the projection of the set X C IP onto the aleatory subspace and 
by Proj e (A) its projection onto the epistemic subspace. Note that V e = Proj e (T?) and 
V a = Proj a ('D). While the aleatory uncertainty is fully prescribed by a PDF whose support 
set is a subset of T> a , the epistemic variable is prescribed by a subset of V e . In this context, 
w(a,e ) = maxj<„{gj(a, e)} is the worst-case requirement function. 

A. Reliability Analysis 

The reliability analysis of a system subject to the design requirements g(a,e) < 0 consists 
of determining 


Range (P[J r (e)]) 


min P[w(-, e ) > 0], maxP[w(-, e) > 0] , 

eeA e eeAe 


(52) 


and calculating the epistemic point (s) that realize this interval’s limits. The range enclosing 
property (23) becomes 

w(a, e, p) < w(a, e) < w(a, e, p). (53) 

As before, p is a partition of the master domain T> = T> a x T> e . A few manipulations lead to 


P [w(a, e, p) > 0] < P[J r (e)\ < P [ w(a , e, p) > 0] , 


( 54 ) 
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which is equivalent to 


1 - ^k(0; e) < P[F(e)\ < 1 - F w { 0; e), (55) 

where F x (-\e) is the CDF of x(- : e ). Lower and upper bounds to F w (-;e), which are 
parametrized by e, are given by 


Fu e) < F w (x ; e) < Fj^x] e). (56) 

Details on the evaluation of Equations (55) and (56) are presented next. Assume that 
p(T>) is a partition of the master domain resulting from the algorithm in Section IV-A. In this 
case the hyper-rectangles comprising the partition are supported in both the epistemic and 
aleatory spaces. Subdivisions in the aleatory dimensions will tighten the probability bounds 
while subdivisions in the epistemic dimensions will tighten the supersets containing the best- 
and worst-case epistemic realizations. The evaluation of the CDF bounds in Equation (56) 
yields 

£ P[Proj„(Kj)]<F„(i ie )< £ P[Proj„(TCj)], (57) 

j£j(x,e) tSJ(x,e) 

where first J(e) = {j : 1 < j < t and e G Proj e (77j)} contains the indices of the members 
of the partition p(T>) whose projection onto the epistemic subspace contains the epistemic 
realization e, and then J(x,e) = {j G J(e) : w(p^) < x} and J_(x,e) = {j G J(e) : 
w(p u) ) < x} pick out those members of the partition which contribute to the respective 
CDFs. Note the similarity between Equations (36) and (57). The evaluation of Equation 
(57) at x — 0 leads to Equation (55). 

The bounds in Equation (57) enable determining an outer bounding interval to the failure 
probability range 


Range (P[J r (e)]) C [1 — maxF„(0; e), 1 — min F w ( 0; e))]. (58) 

e£ A e eeA e 

Since w and W are piecewise constant on a hnite partition of A, the extrema in Equation 
(58) can be determined exactly in a hnite number of steps. 

Supersets of e* and e*, where such extrema occur, are given by 

£* = (J |Proj e (77j) : Fb(0; e) > maxF^O; e)| D e*, (59) 

jeJ( 0,e) 

£* = [J |Proj e (7^j) : FW(0;e) < minF^(0;e)| D e*. (60) 

jel( 0,e) 


The inequalities in these expressions result from combining Equations (27-28) and (55). 
Equation (57) enables evaluating these inequalities. Observe that the subsets Proj e (77j) of 
T> e , can be used to generate a partition of V e into finitely many pieces, on each of which 
Fw( 0; e) and F w ( 0; e) is constant. Equation (58) converges to (52) while £ ifl and £* converge 
to e* and e* when the volume of A approaches zero. 

The algorithmic implementation of the above ideas has the very same structure of the 
algorithm in Section IV-A. Only two minor variations are required. First, at Step 1, 7Z will 
now be the largest element of A* that satisfies Proj e (77) G £* U£*. This will ensure that only 

15 of 25 


American Institute of Aeronautics and Astronautics 



those regions of the epistemic space where the best- and worst-case realizations may be lo- 
cated are partitioned further. Second, at Step 5, the new stopping criterion could be given by 
Volume [£* U £*] being smaller than a prescribed constant, or by having a sufficiently small 
change in the size of the bounding interval in Equation (58). The former criterion implies 
that the supersets are sufficiently small, while the latter implies tight bounds on the extrema. 

Example 3: Consider the requirement functions 

91 = ~3pIpI + 3p? ~p\~ 4 p 2 - 15, (61) 

9 2 = - P l iP-2 + PiPl - ^ (62) 

93 = -2 pi - p\ - p\P2 - (63) 

dehned over the master domain V G [—2, 2] x [—2, 2], 

To start, p l will be considered a Beta-distributed aleatory variable with parameters (4, 4) 
in [—2,2] while p 2 will be epistemic in [—2,2], Therefore, a = {Pi}, e = {p 2 } and A = T>. 
In this setting we would like to bound the extrema of P[J-\ and calculate supersets of the 
corresponding best- and worst- case epistemic values. Figure 4 shows the supersets £* and 
£* of the worst- and best-case epistemic realizations as a function of the iteration number. 
This figure illustrates how the supersets are progressively refined. In both cases, there seems 
to be a single extrema. The existence of an isolated region in the vicinity of p 2 = 0 for £* 
at iteration number 8 and of p 1 = 0.39 for £* at iteration number 14, and their subsequent 
disappearance, suggests the existence of local maxima. Figure 5 shows the corresponding 




Figure 4. 


Sequence of E* (upper graph) and A (lower graph) for a 


{p x } and e = {p 2 }. 


partition of the parameter space. Note that A, the region containing the failure domain 
boundary colored in white, has not been uniformly refined. Further notice that the density 
of boxes increases with the proximity to the extrema of P[J r (e)]. These extrema occur within 
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E* = [1.9336,2] and £* = [—1.07, —0.941] at iteration % = 20. The top of Figure 6 shows the 



Figure 5. Partition of the uncertain parameter space for a = {p x } and e = {p 2 }- 


probability of A as a function of the epistemic variable p 2 ■ The bottom of this figure shows 
the failure probability bounds over the epistemic domain. Note that P[ A] approaches zero 
in the vicinity of the points where the best- and worst- case epistemic values occur. Values 
of p 2 where P[ A] is comparatively large are regions of no interest. The supersets £* and £* 
corresponding to the last iteration in Figure 4 are superimposed in Figure 6. Finally, Figure 




Figure 6. Probability of A and bounds of P[F] for a = {pi} and e = {p 2 }- 


7 shows the ensemble of all CDF bounds, as well as bounds corresponding to epistemic 
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values in £* and £*. While the lower bound of the worst-case representative crosses w = 0 
at Fb(0; e) = 1 — 0.076931, the upper bound of the best-case representative crosses w = 0 
at Fu j(0; e) = 1 — 0.97148. Hence, the bounding interval of the failure probability range 
is [0.076931,0.97148]. Note that the separation between the lower and upper CDF bounds 
decreases with the proximity to w = 0. 



Figure 7. Ensemble of CDF bounds, worst- and best-case CDF bounds representatives for 
a = {Pil and e = {p 2 }- 


We will now repeat this example by considering p x as an epistemic variable while p 2 will be 
an aleatory variable having a Beta PDF with parameters (4, 6) such that A = T>. Therefore 
a = {p 2 } and e = {jq}. Figure 8 shows the supersets £* and £* of the worst- and best-case 
epistemic realizations as a function of the iteration number. In contrast to the previous 
example, the £* sequence suggests the existence of two global minima. Recall that the exact 
nature of the approach guarantees that all global extrema will be found. Figure 9 shows the 
corresponding partition of the parameter space. Note that the density of boxes in A increases 
with the proximity to £* = [1.2578,2] and £* = [-0.48434,-0.46094] U [0.46094,0.48434], 
These are the supersets of the best- and worst-case epistemic realizations at iteration 20. 
Figure 10 and 11 are analogous to Figures 6 and 7. While the lower bound of the worst-case 
representative crosses w = 0 at 0; e) = 1 — 0.01792, the upper bound of the best-case 
representative crosses w = 0 at F, 1; (0; e) = 1 — 0.74623. Hence, the bounding interval of the 
failure probability range is [0.01792,0.74623]. 

B. Performance Analysis 

E[g(a,e)\ and V[g(a,e)\ can be calculated analytically when the aleatory parameters are 
independent random variables and g is a single polynomial. In such a case one can solve 
for the corresponding extrema and the best- and worst- case realizations by performing the 
analysis of Section III to the resulting expressions. These expressions will be polynomials in 
e. If g is a piecewise polynomial function of p (e.g., the worst-case requirements function 
w) and the aleatory parameters are arbitrarily dependent random variables, an approach 
analogous to the one in the preceding section can be applied. 
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Figure 8. Sequence of £* (upper graph) and £* (lower graph) for a = {p 2 } and e = {Pi}- 
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Figure 9. Partition of the uncertain parameter space for a = {p 2 } and e = { p 1 }. 

( 64 ) 

( 65 ) 


An outer bounding interval to the range of expected values is given by 


Range (E[g(-, e)]) C 


min E\g(-, e)], ma: xE[g(-, e)] 

eGAe — eGAe 


where 


%(-,£)]= E s(P U> ) p l Proj^Kj)], 


19 of 25 



American Institute of Aeronautics and Astronautics 






Figure 10. Probability of A and bounds of P[J-\ for a = {p 2 } and e = {p^- 



Figure 11. Ensemble of CDF bounds, worst- and best-case CDF bounds representatives for 
a = {P-i} and e = {pj. 


E[g(;e)]= J2s(P U) ) p l Pr°ia(^)], (66) 

j6J(e) 

and J(e) = {j : 1 < j < t, e G Proj e (7?.j)}. Supersets of the epistemic realizations prescrib- 
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ing the endpoints of E[g(-, e)] are 


£*= IJ { Pr °je(^) : -%(-> e )] < mm-%]} D e*, 

ieJ(e) 

(67) 

£* = IJ { Pr °je(^) : E [d(-,e)] > max^^j D e* 

ieJ(e) 

(68) 

Equations (65) and (66) enable evaluating the inequality constraints in 
(68). 

An outer bounding interval of the variance is 

Equations (67) and 

Range (V\g(-, e)]) C [max{0, a(e)}, /3(e)] , 

(69) 


where a(e) and /3(e) are given by Equations (50) and (51) with the expected values given 
by Equations (65) and (66), and with /i(e) and u(e) given by 

M<0= V g(p u) ) 2 P [Proj 0 (TC,)] + ^ J(p«>) 2 P[ Proj„(TC,)], (70) 

jeU(e) jev (e) 

K e ) = ^(p 0) ) 2p [ pr °ja(^)] + g(p U) ) 2 P[Pro ) a (TZj)}, (71) 

j£U(e)UW(e) jeV(e) UX(e) 

where U, V, W and A" are defined in Equations (45-48) and e e Proj e (77j). Supersets of 
the epistcmic realizations where the extreme values of V[g(-,e)\ occur are 

£* = |Proj e (77j) : max{0, ce(e)} < min/3(e) j D e *, (72) 

ieJ(e) 

£* = |Proj e (77j) : /3(e) > max {max{0, a(e)}}| D e*. (73) 

ieJ(e) 

These expressions result from applying Equations (27-28) to the bounding functions in (69). 
As before, even though the approach is applicable to any partition of the master domain, 
tighter bounds are obtained if the partition is hirer where where either g — g or E[g — g\Tl] 
are large. 


VI. Global Sensitivity Analysis 

In the context of UQ, sampling-based methods 4 and analysis of variance are widely 
used sensitivity analysis techniques. These methods evaluate sensitivities by quantifying the 
contribution of the uncertainty in a particular parameter to the spread in the variance of a 
performance function. While this information provides some insight as to what parameters 
are dominant and which ones are not, it fails to account for the effect of such uncertainty in 
other statistics. Therefore, it is possible that the contribution of the uncertainty in a given 
parameter to the variance is insignificant, while its contribution to the failure probability is 
substantial. 

In this section we consider the sensitivity of the mean, variance and failure probability to 
the uncertainty model of p. The sensitivities are global since they are not restricted to a single 

21 of 25 


American Institute of Aeronautics and Astronautics 



point of the uncertain parameter space, nor to a particular family of distribution functions. 
If the sensitivity is small, the statistic will not change noticeably when the uncertainty model 
is changed within the support set, e.g., the UQ of a system where an unimportant parameter 
is modeled as a random variable supported in A is practically identical to the UQ of the 
system where such a parameter takes any fixed value in A. This practice, which can be 
used to reduce the number of parameters that should be considered uncertain, facilitates the 
deployment and convergence of UQ methods whose performance, accuracy and correctness is 
critically linked to the dimension of the parameter space, e.g., polynomial chaos, homothetic 
deformations, FORM, etc. 

For instance, consider the reliability analysis of a system that depends on s independently 
distributed aleatory variables. Our objective is to determine the sensitivity of the failure 
probability P[w{jp) > 0] to the uncertainty model of p l . In principle, for each of the infinite 
number of models that could be chosen to describe p x there is a different value of the failure 
probability. Note that failure probability range that results from propagating all possible 
uncertainty models of p 1 is bounded by the failure probability range corresponding to the 
subfamily of Dirac deltas. Consequently, we will describe the aleatory variable p x as if it 
were epistcmic. In this setting, the methods of Section V can be used to calculate an outer 
bounding interval to the failure probability range for a = {p 2 , . . . ,p s } and e = {Pi} e A e . 
Note that the failure probability corresponding to any uncertainty model of p 1 supported in 
A e , -even a probabilistic one- will fall within this interval. Therefore, if the bounding interval 
is sufficiently small, all uncertainty models of p 1 supported in A e will lead to practically the 
same failure probability value. By assuming that p l takes on a fixed, known value in A e 
without changing the model of \p 2 , . . . ,p s \, we will obtain a failure probability value that 
differs from any other one we could have obtained by less than the spread of the bounding 
interval. Therefore, small bounding intervals will indicate the p 1 can be safely assumed to 
take on a constant value. If the spread of the bounding interval is large and the partition 
of A near e* and e* is sufficiently fine, p l is a dominant parameter. In this case the failure 
probability will be highly sensitive to the uncertainty model of p 1 . 

A mathematical framework supporting these ideas is presented next. Without loss of 
generality, consider the statistics E[w\, V[w\ and R = P[w > 0] = P[P}- Assume that lower 
and upper bounds to this function, namely w(p) and w(p), are calculated via Bernstein 
polynomials. Let L be a set with all the indices of the parameters in p whose sensitivity is 
to be determined, while K contains the indices of the remaining aleatory variables. Further 
assume that the uncertainty model of the variables in L is independent of the uncertainty 
model of the variables in K. L, which can contain indices of both epistcmic and aleatory 
variables, has m elements where 0 < m < p — 1. Denote by A q the support set of q — 
[p Ll , ■ ■ . , P/\ r J, and by A. r the support set of r e R s ~ m . 


Our objective is to calculate the sensitivity of the statistics in 

Equations (6-8), given by 

E [w; L] (q) = j 

f w(q,r)f r (r)dr, 
A r 

(74) 

R [w, L] (. q ) = j 

( f r (r )dr , 

w(q,r)>0 

(75) 

V [w; L\ {q) = E[w 2 -L]-E[w-L}\ 

(76) 
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to the uncertainty model of the variables in q. In this setting, the global sensitivities are: 

a R (w,w, L) = rna xE[ui] L] — min E[w ; L\, (77) 

qG Aq qG A, 

ctr(w , w, L ) = max R\w] L] — min R[w: L] , (78) 

q£A g qGAg 

ay(w, w, L) = ma xu[ur, ui] L] — min l[w] ui] L ], (79) 

qG Aq qG Aq 

where u and l are given by the interval equation 

| l[w; w ; L], u[w; fa; L]] = E [[w, w] 2 ; L] — L], E[iu; L]] 2 , (80) 

set according to the developments in Section IV-B. These expressions can be calculated 
using the developments of Section V by treating q as if it were an epistcmic variable. 

If a global sensitivity is e, the corresponding statistic will not change more than e when 
we chose any uncertainty model for q supported in A q . When this model is comprised by a 
collection of Dirac deltas, i.e., when the parameters take on fixed, known values; an accurate 
reduction in the uncertain parameter dimension is attained. In general, the sensitivities 
aE(w,w,L) } <jr(w,w, L), and ay(w, w, L) are, in general, independent of each other; e.g., 
some parameters may be unimportant according to the failure probability but they may be 
dominant according to the mean. These following example illustrates these ideas quantita- 
tively. 


Example 4: Consider the requirement functions 


9 1 

92 


4 p\ -p'i + Pi 

1 3 2 1 

-PlP? + “ 

5 12 1000 


- 4 pi 

p\pI- 


400 p\ 


400 pi 


1 S 


PiPl 


m Pl + 2 Pl 2’ 


(81) 

(82) 


defined over the master domain T> G [—2,2] x [—2,2], Assume that p l is a Beta-distributed 
random variable with parameters (7,7), and p 2 is a Beta-distributed random variable with 
parameters (1, 10) both supported in [—2, 2], We want to determine the global sensitivity of 
E[w], V[w\ and P[w > 0] = P[g > 0 ], where w is the worst-case requirement function, to 
the uncertainty models of p 1 and p 2 . 

While the sensitivities for the setting q = {p x } are a E (w,w,l) = —2.20391— (—2.79632) = 
0.5924, a v (w,w, 1) = 1.09538 - 0 = 1.09538, and a R (w,w, 1) = 0.0042194 - 0 = 0.0042194; 
those for q = {p 2 } are <te(w,w, 2) = —0.56841 — (—2.57057) = 2.0022, ay(w,w,2) = 
2.86507 - 0 = 2.86507, and a R (w,w, 2) = 0.47571 - 0 = 0.47571. Figure 12 shows the 
bounds on the mean and variance functions from where the sensitivities were calculated. 
The largest and smallest values of the bounding functions determine o R and ay. While the 
sensitivity of the failure probability to uncertainty in p 2 is 112 times larger than that for p { , 
the sensitivity of the mean to uncertainty in p 2 is 3.38 times larger than that for p 1 . If a 
failure probability of less than 0.0042194 is acceptable, p x can be assumed to take on any 
constant value in [—2,2], 

In this example, estimates of the functions E[w(q)\ and V[w(q)\ can be readily calculated 
using Monte Carlo sampling. These functions will fall between the bounds shown in Figure 
12. The extrema of such functions, which can be calculated via optimization, can be used 
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to approximate the global sensitivities above. This practice will render meaningful values as 
long as the optimizer has converged to the global extrema. Unfortunately, the moments of 
piecewise polynomial functions may not only have derivative discontinuities but also multiple 
local extrema. In contrast to this practice, the approach proposed yields results whose 
correctness is formally verifiable. 






q=P, q=p 2 


Figure 12. Mean (top) and variance (bottom) bounds for q = {jq} (left) and q = {p 2 } (right). 


VII. Conclusions 

This paper presents an uncertainty quantification framework that enables the rigorous 
analysis of polynomial systems subject to both aleatory and epistemic uncertainties. The 
Bernstein expansion approach enables the calculation of analytical bounds to moments and 
failures probabilities as well as finding supersets of the corresponding best- and worst-case 
epistemic realizations. These bounds and supersets can be made as tight and as small 
as desired. Furthermore, the framework enables determining the importance of particular 
uncertain parameters according to the manner in which they affect the failure probability, the 
mean and the variance of a response metric. This information enables ranking the uncertain 
parameters according to the sensitivity of these three statistics. By modelling inconsequential 
parameters as deterministic quantities, the size of the uncertain parameter space can be 
reduced without incurring in significant error. The method proposed generates bounds to this 
error. The analytical nature of the approach eliminates the numerical error that characterizes 
the sampling-based techniques commonly used to propagate aleatory uncertainties as well 
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as the possibility of under predicting the range of the statistic of interest that may result 
from searching for the best- and worst-case epistemic values via nonlinear optimization or 
sampling. 
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